Size Spectrum Species Sensitivity

Exploring Spectrum Sensitivies to Species Omission

For the sensitivity testing we estimated the size spectrum slope and intercept for every year and within each survey area an additional time, at each pass one of the species was omitted from the analysis. This was done to capture the relative importance each species had on the estimation of the biomass size spectra.

Loading Bio Data

The biological data used as a starting point is the same “node” used for the size spectrum estimations using the full community. This is loaded so that we can interrogate the overall biomass patterns that correspond to the results we see from the sensitivity testing.

# Access Biological Data with {gmRi}
# nmfs_bio <- gmRi::gmri_survdat_prep(survdat_source = "bio")
withr::with_dir(rprojroot::find_root('_targets.R'), 
                tar_load(nefsc_1g_binned))


# Get Yearly and Regionally averaged Biomasses
yr_szn_overall <- nefsc_1g_binned %>% 
  group_by(Year, survey_area) %>% 
  summarise(overall_bio = sum(strat_total_lwbio_s, na.rm = T),
            overall_abund = sum(strat_total_abund_s, na.rm = T),
            .groups = "drop")

# Get what it is for each species
sp_yr_szn <- nefsc_1g_binned %>% 
  group_by(comname, spec_class, fishery, hare_group, Year, survey_area) %>% 
  summarise(sp_bio = sum(strat_total_lwbio_s, na.rm = T),
            sp_abund = sum(strat_total_abund_s, na.rm = T),
            .groups = "drop")


# Join in the overall, get the fraction and what the percent would be
perc_bio <- left_join(sp_yr_szn, yr_szn_overall, by = c("Year", "survey_area")) %>% 
  mutate(
    bio_frac = sp_bio / overall_bio,
    abund_frac = sp_abund / overall_abund,
    bio_perc = bio_frac * 100,
    abund_perc = abund_frac * 100
  )

Overall Biomass Distributions

# perc_bio overall rankings
region_overalls <- perc_bio %>% 
  group_by(comname, spec_class, fishery, hare_group, survey_area) %>% 
  summarise(avg_bio_perc   = mean(bio_perc, na.rm = T),
            avg_abund_perc = mean(abund_perc, na.rm = T),
            .groups = "drop")

Biomass

# Treemaps
ggplot(region_overalls) +
  geom_treemap(aes(area = avg_bio_perc, fill = comname)) +
  geom_treemap_text(aes(area = avg_bio_perc, fill = comname, label = comname)) +
  facet_wrap(~survey_area, ncol = 2) +
  scale_fill_gmri() + 
  theme(legend.position = "none") +
  labs(title = "Fraction of Stratified Biomass: 1970-2019")

Abundance

# Treemaps
ggplot(region_overalls) +
  geom_treemap(aes(area = avg_abund_perc, fill = comname)) +
  geom_treemap_text(aes(area = avg_abund_perc, fill = comname, label = comname)) +
  facet_wrap(~survey_area, ncol = 2) +
  scale_fill_gmri() + 
  theme(legend.position = "none") +
  labs(title = "Fraction of Stratified Abundance: 1970-2019")

Yearly Patterns

Biomass

# Play around with plotting yearly influences
ggplot(perc_bio, aes(Year, bio_perc)) +
  # geom_point(aes(color = comname), show.legend = F) +
  # geom_line(aes(color = comname, group = comname), show.legend = F) +
  geom_col(aes(fill = spec_class)) +
  facet_grid(survey_area~.) +
  scale_color_gmri() +
  scale_fill_gmri() +
  theme(legend.position = "bottom") +
  labs(y = "Percent of Yearly Biomass",
       x = NULL,
       fill = "")

# 

Abundance

# Play around with plotting yearly influences
ggplot(perc_bio, aes(Year, abund_perc)) +
  # geom_point(aes(color = comname), show.legend = F) +
  # geom_line(aes(color = comname, group = comname), show.legend = F) +
  geom_col(aes(fill = spec_class)) +
  facet_grid(survey_area~.) +
  #facet_grid(spec_class~survey_area) +
  scale_color_gmri() +
  scale_fill_gmri() +
  theme(legend.position = "bottom") +
  labs(y = "Percent of Yearly Biomass",
       x = NULL,
       fill = "")

# 

Loading Sensitivity Results

# Access Biological Data with {gmRi}
# nmfs_bio <- gmRi::gmri_survdat_prep(survdat_source = "bio")
withr::with_dir(rprojroot::find_root('_targets.R'), 
                tar_load(species_sensitivity_shifts))

# Make year a numeric value
sens_results <- species_sensitivity_shifts %>%
  mutate(Year = as.numeric(Year)) 


# Join in the biomass info
sens_results <- sens_results %>% 
  select(-spec_class) %>% 
  left_join(perc_bio,
            by = c("Year" = "Year", 
                   "survey_area" = "survey_area", 
                   "spec_omit" = "comname"))

Single Species Removal Impacts

sens_ranks <- sens_results %>% 
  group_by(
    `Species Removed` = spec_omit, 
    Area = survey_area, 
    Guild = spec_class, 
    Fishery = fishery) %>% 
  summarise(`Avg. Slope Shift` = mean(slope_shift, na.rm = T),
            `Avg. Slope Shift (z)` = mean(slope_shift_z, na.rm = T),
            .groups = "drop") %>% 
  mutate(across(where(is.numeric), round, 2))

Slope Increase

# Slope Steepeners
sens_ranks %>% 
  arrange(`Avg. Slope Shift (z)`) %>% 
  head(8) %>% 
  gt::gt()  %>%
  gt::tab_header(
    title = "Top 8 Steepeners",
    subtitle = "Removing these species causes slope to steepen"
  )

Top 8 Steepeners
Removing these species causes slope to steepen
Species Removed	Area	Guild	Fishery	Avg. Slope Shift	Avg. Slope Shift (z)
bullnose ray	SNE	NA	Not Labelled	-0.41	-2.16
silver hake	GoM	Groundfish	Commercially Targeted	-0.20	-1.08
silver hake	SNE	Groundfish	Commercially Targeted	-0.14	-0.75
barndoor skate	GB	Elasmobranch	Commercially Targeted	-0.11	-0.59
butterfish	MAB	Pelagic	Not Commercially Targeted	-0.08	-0.58
silver hake	GB	Groundfish	Commercially Targeted	-0.11	-0.56
roughtail stingray	SNE	NA	Not Labelled	-0.09	-0.49
spotted hake	MAB	Groundfish	Not Commercially Targeted	-0.06	-0.44

Slope Flattener

# Slope Flatteners
sens_ranks %>% 
  arrange(desc(`Avg. Slope Shift (z)`)) %>% 
  head(8) %>% 
  gt::gt()  %>%
  gt::tab_header(
    title = "Top 8 Flatteners",
    subtitle = "Removing these species causes slope to flatten"
  )

Top 8 Flatteners
Removing these species causes slope to flatten
Species Removed	Area	Guild	Fishery	Avg. Slope Shift	Avg. Slope Shift (z)
spiny dogfish	MAB	Elasmobranch	Commercially Targeted	0.13	0.92
spiny dogfish	SNE	Elasmobranch	Commercially Targeted	0.15	0.80
winter skate	GB	Elasmobranch	Commercially Targeted	0.04	0.22
spiny dogfish	GoM	Elasmobranch	Commercially Targeted	0.04	0.21
atlantic cod	GB	Groundfish	Commercially Targeted	0.04	0.19
spiny dogfish	GB	Elasmobranch	Commercially Targeted	0.03	0.17
atlantic cod	GoM	Groundfish	Commercially Targeted	0.02	0.12
sand tiger	MAB	NA	Not Labelled	0.02	0.12

Regional Influence Heatmaps

# Do we need to show the relative influence?

# Things seem to wash out with the raw z-scores
# Maybe they should be adjusted/weighted for how much any single omission
# pulls a given year compared to either the 
# average z-shift for that year+region
# or by the max that it shifts....

# plot heatmap
influence_heatmap <- function(sens_results, 
                              filt_area, 
                              lims = c(-1, 1),
                              rescale_z = F){
  
  # Make breaks to indicate limits are truncating
  blabs <- seq(from = lims[1], to = lims[2], length.out = 5)
  blabs[c(1,5)] <- c(str_c(">", blabs[1]), 
                     str_c("<", blabs[5]))
  
  # Filter data
  filt_dat <- sens_results %>%
    mutate(Year = as.numeric(Year)) %>%
    filter(survey_area == filt_area)
  
  # Rescale?
  if(rescale_z == T){
    message("Rescaling not finished")
    
    
  }
                     
  
  # Make heatmap
  filt_dat %>%
  ggplot(aes(Year, fct_rev(spec_omit), fill = slope_shift_z)) +
    geom_tile() +
    scale_fill_distiller(palette = "RdBu",
                         limits = lims, 
                         labels = blabs,
                         oob = scales::oob_squish, 
                         na.value = "gray20") +
    scale_x_continuous(expand = expansion(add = c(0,0))) +
    facet_wrap(~survey_area, ncol = 2, scales = "free") +
    #facet_grid(spec_class~survey_area, scales = "free") +
    #facet_grid(fishery~survey_area, scales = "free") +
    theme(legend.position = "bottom", 
          legend.title = element_text(angle = 0), 
          axis.text.y = element_text(size = 7),
          panel.background = element_rect(fill = "gray40"),
          panel.grid = element_line(color = "transparent")) +
    guides(fill = guide_colorbar(title.position = "top", 
                                 title.hjust = 0.5, 
                                 barwidth = unit(3.25, "in"))) +
    labs(y = "Species omitted", 
         fill = "Shift in Slope (z-score)")
  
  
}

Gulf of Maine

# GOM
sens_results %>% 
  #filter(spec_class == "Groundfish") %>% 
influence_heatmap(filt_area =  "GoM", lims = c(-2,2))

Georges Bank

# GOM
influence_heatmap(sens_results, "GB", lims = c(-2,2))

Southern New England

# GOM
influence_heatmap(sens_results, "SNE", lims = c(-4,4))

Mid-Atlantic Bight

# GOM
influence_heatmap(sens_results, "MAB", lims = c(-4,4))

Correlation to Biomass Fraction

sens_results %>% 
  ggplot(aes(bio_frac, slope_shift_z)) +
  geom_point(alpha = 0.3) +
  geom_smooth(formula = "y ~ x", method = "lm") +
  labs(x = "Fration of Biomass", y = "Shift in Spectra Slope (z)")

Regional Influence Horizon Plots

I think this might help get around the “wash out” of the heatmap. The color scale should be longer, or potentially just not center on white…

Singe Species Checks:

# plot a species:
# multiple panel plot function
# Function for plotting a species that stands out
plot_spectra_influence <- function(sens_results, filter_spec, filter_area){

  # Grab the relevant data
  filt_d <- sens_results %>%
    filter(spec_omit == filter_spec,
           survey_area == filter_area)
  
  # Make Plot Label
  p_label <- str_c(str_to_title(filter_spec), " Omission")
  
  
  # 1. Actual Value Shift
  # 1a. Slope
  splot <- filt_d %>%
    ggplot(aes(Year, l10_slope_strat)) +
    geom_line(aes(color = "All Species"), group = 1) +
    geom_point(aes(color = "All Species")) +
    geom_line(aes(y = omit_slope_strat, color = p_label), group = 1) +
    geom_point(aes(y = omit_slope_strat, color = p_label)) +
    facet_wrap(~survey_area, ncol = 1) +
    scale_color_gmri() +
    theme(legend.position = "bottom") +
    labs(color = NULL,
         y = "Biomass Spectrum Slope",
         x = NULL)
  
  # 2. What the z-score shift is
  z_shift <- filt_d %>% 
    ggplot(aes(Year, slope_shift_z)) +
    geom_line() +
    geom_point() +
    labs(y = "Slope Shift (z)",
         x = NULL)
  
  # Biomass Fraction Each Year:
  bfrac <- filt_d %>%
    ggplot(aes(Year, bio_frac)) +
    geom_line() +
    geom_point() +
    scale_y_continuous(labels = scales::percent) +
    labs(y = "Biomass Fraction")
  
  
  # Do they fit
  plotsemble <- splot / z_shift / bfrac
  return(plotsemble)
}

Haddock (nowhere)

plot_spectra_influence(sens_results = sens_results, 
                       filter_spec = "haddock", 
                       filter_area = "GoM")

Spiny Dogfish

plot_spectra_influence(sens_results = sens_results, 
                       filter_spec = "atlantic herring", 
                       filter_area = "SNE")

Silver Hake (GoM & GB)

plot_spectra_influence(sens_results = sens_results, 
                       filter_spec = "silver hake", 
                       filter_area = "GoM")

plot_spectra_influence(sens_results = sens_results, 
                       filter_spec = "silver hake", 
                       filter_area = "GB")

Butterfish

plot_spectra_influence(sens_results = sens_results, 
                       filter_spec = "butterfish", 
                       filter_area = "GoM")

Atlantic Cod

plot_spectra_influence(sens_results = sens_results, 
                       filter_spec = "atlantic cod", 
                       filter_area = "GoM")